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Abstract 

This paper establishes a remarkable result regarding Palm distribu¬ 
tions for a log Gaussian Cox process: the reduced Palm distribution 
for a log Gaussian Cox process is itself a log Gaussian Cox process 
which only differs from the original log Gaussian Cox process in the 
intensity function. This new result is used to study functional sum¬ 
maries for log Gaussian Cox processes. 

Keywords: J-function; joint intensities; Laplace approximation; nearest- 
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1 Introduction 


Palm distributions (see e.g. Mpllcr and Waagepetersen, 2004 Daley and Vere- 


Jones, 2008) are important in the theory and application of spatial point 


processes. Intuitively speaking, for a prespecified location in space, the Palm 
distribution of a point process, with respect to this location, plays the role 
of the conditional distribution of the point process given that the aforemen¬ 
tioned location is occupied by a point of the point process. 
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The present paper focuses on log Gaussian Cox processes (Mpllcr et al. 


1998) which provide very flexible, useful, and popular models for modeling 
spatial patterns in e.g. biology and spatial epidemiology. The paper estab¬ 
lishes a surprisingly simple characterization of Palm distributions for such a 
process: The reduced n-point Palm distribution is for any n > 1 itself a log 
Gaussian Cox process that only differs from the original log Gaussian Cox 
process in its intensity function (not to be confused with the random intensity 
function generating this kind of Cox process). This result can be exploited 
for functional summaries as discussed later. The simplicity and completeness 
of this result is remarkable when compared with Palm distributions for other 
common classes of spatial point processes. Reduced Palm distributions for 
Gibbs point processes are also themselves Gibbs point processes but with 
densities only known up to a normalizing constant. For shot-noise Cox pro¬ 


cesses (Mpllcr, 2003) one-point reduced Palm distributions have a simple 


characterization as cluster processes similar to shot-noise Cox processes but 
this is not the case for n-point Palm distributions when n > 1. 

The paper is organized as follows. Section [2] reviews the general definition 
of reduced Palm distributions of any order and relates this to Cox processes. 
Section [3] establishes our characterization result for log Gaussian Cox pro¬ 
cesses. Section [4] applies this result to functional summaries for stationary log 
Gaussian Cox processes, in particular the so-called F, G, and J-functions, 
where we establish some new theoretical results, consider how to calculate 
F, G , and J using Laplace approximations, and discuss an application. Sec¬ 
tion [5] concludes the paper. 

2 Palm distributions 


Our general setting is as follows. For ease of exposition we view a point 
process as a random locally finite subset X of a Borel set S C M d , d > 1; 


for measure theoretical details, see e.g. Mpller and Waagepetersen (2004) or 


Daley and Vere-Jones (2003). Denoting = X n B the restriction of X 


to a set B C S, local finiteness of X means that X# is finite almost surely 
(a.s.) whenever B is bounded. We denote J\f the state space consisting of 
the locally finite subsets (or point configurations) of S. We use the generic 
notation h for an arbitrary non-negative measurable function defined on A/”, 
S n , or S n x J\T for n = 1, 2,.... Furthermore, Bo is the family of all bounded 
Borel subsets of S. Finally, recall that the void probabilities P(X# = 0), 
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K C S compact, uniquely determine the distribution of X. 


2.1 Factorial moment measures and Palm distributions 


This section provides the general definition of reduced Palm distributions of 
any order. For finite point processes specified by a density, a simpler and 
more explicit definition is available as reviewed in Coeurjolly et al. (2015). 

For n = 1,2,... and Bi e Bo, the n-th order factorial moment measure 
a( n ' > is defined by 


(B\ x i ?2 x ■ ■ • x B n ) = E 1 [x\ e B \,..., x n 6 B n ), 

x\ 


where l(-) denotes the indicator function and 7 ^ over the summation sign 
means that X\,... ,x n are pairwise distinct. If a ^ has a density p ^ with 
respect to Lebesgue measure, p^ n > is called the n-th order joint intensity 
function and is determined up to a Lebesgue nullset. Therefore, we can as¬ 
sume that p( n \x 1 ,..., x n ) is invariant under permutations of X\, ..., x n , and 
we need only to consider the case where x \,..., x n G S are pairwise distinct. 
Then p^ n \x 1 ,..., x n ) d^i • ■ • dx n can be interpreted as the approximate prob¬ 
ability for X having a point in each of infinitesimally small regions around 
xi,... ,x n of volumes day,... dx n , respectively. We also write p{u) for the 
intensity function pd)(w). 

Moreover, for any measurable F C J\f , define the n-th order reduced 
Campbell measure C^ n >' as the measure on S n x J\f given by 

C {n) \B x x B 2 x ••• x B n x F) = 

E ^2 1(^1 e B u ... ,x n e B n ,X \ {xi,... ,x n } e F). 


Note that x F), as a measure on S n , is absolutely continuous with 

respect to a^ n \ with a density P' £i (F) which is determined up to an a^ 
nullset, and a^HAf) = C ( - n ^{B\ x B 2 x ■ ■ ■ x B n x A f). By the so-called 
Campbell-Mecke formula/theorem, we can assume that P' Xl Xn ( m ) a P°i n f 
process distribution on J\f, called the n-th order reduced Palm distribution 

Daley and Vere-Jones, 2008). We denote by X^ x 

Again we need only to 


given xi,... ,x n (see e.g. 


a point process distributed according to P x 


i — 1 - 

Xi,...,X n ' 
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consider the case where xi,... ,x n are pairwise distinct. Then P' xl ,..., Xn can 
be interpreted as the conditional distribution of X \ {xi, given that 

Xi,..., x n G X 

If phO exists, then by standard measure theoretical arguments we obtain 
the extended Campbell-Mecke formula 




E ^ h(x i,...,x„,X\{xi,...,x n }) 

= / ■■■ / Eh(xi,... ^X^ ,„ !Xn )p (n) (x 1 ,... ,x n ) dxi • ■ • dx n . (1) 


Suppose p( m+? d exists for an m > 1 and n > 1. Then, for pairwise distinct 
«i,..., u m , xi ,..., x n G S, it follows easily by expressing a^ m+ri ^ as an expec¬ 
tation of the form ({TJ) that X; E] Xn has m-th order joint intensity function 


p (m) i 
rxi,...,x n 


Ui, 






0 


if p (n) (xi,. ..,x n ) > 0, 
otherwise. 


( 2 ) 


We also write p xl ..., Xn f° r the intensity function pxl...x ri - 


2.2 Cox processes 


Let A = {A(x)} a . eS be a nonnegative random held such that A is locally 
integrable a.s., that is, for any B e B 0 , the integral f B A(x)dx exists and 
is finite a.s. In the sequel, X conditional on A is assumed to be a Poisson 
process with intensity function A; we say that X is a Cox process driven 
by A. We also assume that A has moments of any order n — 1, 2,.... Then 
the joint intensities of X exist: For any n = 1,2,... and pairwise distinct 
X\i • • •, x n G S , 


P {n) {x 1 ,.. ■, x n ) = E 


n A( ^) 


(3) 


The following lemma, which is verified in Appendix A, gives a character¬ 
ization of the reduced Palm distributions and their void probabilities. 
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Lemma 1. Let X be a Cox process satisfying the conditions above. Then, 
for any n = 1,2,..., pairwise distinct x\,,x n £ S, and compact K C S, 


p {n \xi ,..., x n )E {h (a*, ...,x n , } = E \ K x u ■ ■ ■, x n , X) Y[ A (xi) 


2—1 


and 

p {n \x±,... ,x n )v{yf xi ^ Xn n/i = 0) = E 


(4) 


exp|— J A(u) drt j> riA(^) 


(5) 


3 Reduced Palm distributions for log Gaus¬ 
sian Cox processes 

For the remainder of this paper, let X be a Cox process driven by A = 
{A(x)} xS 5 , where A(x) = exp{Y(x)} and Y = {Y(x )} a;e 5 is a Gaussian 
random field with mean function p and covariance function c so that A is 


locally integrable a.s. (simple conditions ensuring this are given in Mpller 


et al. 1998). Then X is a log Gaussian Cox process (LGCP) as introduced 


by Coles and Jones (1991) in astronomy and independently by Mpller et ah 
(1998) in statistics. 

For distinct x,y G S, define the so-called pair correlation function g(x, y) = 
p ( ' 2 \x,y)/{p(x)p(y)} (the following result shows that p > 0 in the present 
case). By Mpller et ah (1998, Theorem 1), 

p(x) = exp{p(x)+c(x,x)/2}, g(x,y) = ex.p{c(x,y)}, x,yeS,x^y, 

( 6 ) 

and for pairwise distinct Xi ,..., x n G S, 


p^ n \x 1 ,...,x n ) = <{ JJp(xi) Yl d( x i, x j) 

„ Kkj<n 


. i=l 


is strictly positive. 

For u,xi,... ,x n E S, define 


Px u ...,xM = p{u) + Y c{u, x d- 


2—1 


(7) 
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Combining (J2| and ©-0. we obtain for any pairwise distinct 
Mi,, u m , x\,..., x n G S with m > 0 and n > 0, 


P. 


M 


(Ml, 


! ^m. ) 


Px 1 ,..., 


1=1 



JJ g(ui,uj)\, (8) 

l<i<j<m 


where 

Px = exp {/i Xl) ...,*„(«) + c(m,m)/2}. 
Thereby the following proposition follows. 


Proposition 1. For the LGCP X and any pairwise distinct xi,...,x n G 
S, X!,. Xn has m-th order joint intensity ([8]) which agrees with the m-th 
order joint intensity function for an LGCP with mean function p xu ... jXn an d 
covariance function c for the underlying Gaussian random field. 


Proposition [l indicates that also XL Xn conlcl be an LGCP. A sufficient 
condition, considered by Macchi (1975), is the existence of a number a = 
a(B ) > 1 for each set B e B 0 such that 


E 


exp \a A (it) d u 


' B 


< oo. 


(9) 


However, we have not been successful in verifying this condition which seems 
too strong to hold for any of the covariance function models we have con¬ 
sidered, including when c is constant (then X is a mixed Poisson process) 
or weaker cases of correlation, e.g. if c is a stationary exponential covariance 
function. The case where c is constant is closely related to the log normal 
distribution which is not uniquely determined by its moments (Heyde, 1963). 

Accordingly we use instead Lemma [l] when establishing the following 
theorem, which implies that the LGCPs X and X^ Xn share the same pair 
correlation function and differ only in their intensity functions. 


Theorem 1. For pairwise distinct aq,..., x n G S, X!,. n is an LGCP with 
underlying Gaussian random field Y x ^^ Xn , whereY xl __ Xn has mean function 
Pxi,...,x„ an d covariance function c. 

Let Y = Y — fi be the centered Gaussian random field with covariance 
function c. Theorem [l] is a consequence of the fact that the probability 
measure of Y Ilvi . A is absolutely continuous with respect to the one of Y, 
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with density exp |X )" =1 v( x i) ~ Y^ij=\ c ( x i, x j)/ 2 j when y is a realization of 

Y. This result is related to the Cameron-Martin-Girsanov formula for one¬ 
dimensional Gaussian processes. A short selfcontained proof covering our 
spatial setting is given in Appendix A. 

Often we consider a non-negative covariance function c or equivalently 
9> 1, which is interpreted as ‘attractiveness of the LGCP at all ranges’, but 
even more can be said: A coupling between X and Xn is obtained by 

taking Y xl> „. )Xn (x) = Y(x) + YTi=i c ( x , x d- Thus, if c > 0 and we are given 
pairwise distinct points x \,..., x n G S, we can consider X as being included 
in X^ Xn , since X can be obtained by an independent thinning of X^ Xn , 
with inclusion probabilities exp{— ^" =1 c(x, £*)}, x G X \ {aq,..., x n }. This 
property clearly shows the attractiveness of the LGCP if c > 0 (equivalently 
9>l). 


4 Functional summaries for stationary log Gaus¬ 
sian Cox processes 


Throughout this section, let S = R d and assume that the LGCP X is sta¬ 
tionary, i.e., its distribution is invariant under translations in W 1 . By (J 6 ])- 
([7]), this is equivalent to stationarity of the underlying Gaussian random 
field Y, that is, the intensity p is constant and the pair correlation function 
g(x,y ) = g(x — y ) is translation invariant, where x, y G are distinct, and 
g(x) = exp{c(a;)} and c(x) = c(o,x ) for x G M d , where o denotes the origin 
in It is custom to call P Q the reduced Palm distribution at a typical 
point, noticing that for any x G R d , X(, and Xj,, — x = {y — x \ y E X^,} are 
identically distributed. 

Denote B (o, r) the ball in of radius r > 0 and centered at o. Popular 
tools for exploratory purposes as well as model fitting and model checking 
are based on the following functional summaries where r > 0 (see e.g 
and Waagepetersen, 2004): 


Mpller 


(i) the pair correlation function g and the related Ripley’s A'-function 
given by 

A» = -E#{X' nB(o,r)}= I g(x). 

P J B(o,r) 

Thus, pK(r) is the expected number of further points in X within 
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distance r of a typical point in X. If g{x) depends only on the distance 
||x|| then g and K are in one-to-one correspondence; 

(ii) the empty space function given by 

F(r) = P{Xn B{o,r) ± 0}, 

which is the probability that X has a point within distance r of an 
arbitrary fixed location; 

(iii) the nearest-neighbour distribution function given by 

G(r) = P{X|,nfl(o,r)/0}, 

which is the probability that X has a further point within distance r 
of a typical point in X; 

(iv) the J-function given by 


J(r) 


1 - G(r) 

1 -F(r)’ 


with the convention a/0 = 0 for any a > 0. 


Section 14.11 establishes some new results for these theoretical functions 
and Section 42 discusses how they can be calculated using a Laplace ap¬ 
proximation. Section 43 illustrates this calculation and Section 44 discusses 
an application for a real dataset. 


4.1 New formulae for G and J 

By conditioning on Y, we see that 


1 — F(r) = E ( exp 


' B(o,r) 


exp{Y'(x)} dx 


( 10 ) 


Using the Slivnyak-Mecke formula, Mpllcr et al. (1998) showed that 


1 - G(r) = - E exp 
P 


Y(o) — / exp{Y(a:)} Ax 

JB{o,r) 


( 11 ) 















Since the nearest-neighbour distribution function for X is the same as the 
empty space function for X;,, which is an LGCP with underlying Gaussian 
random held Y 0 (x) = Y(x) + c(x), and since g(x) = exp{c(a;)}, we obtain an 
alternative expression 


1 - G(r) = E exp 


/ g(x) exp{y(a:)} dx 
J B(o,r) 

Therefore, we also obtain a new expression for the J-function, 


J(r) = 


( 12 ) 


E ^exp 

~ I B (o,r)9(x)exp{Y(x)}dx 

]) 

E 

exp 

— fB(o,r) ex p{V()r)} dx 

) 



(13) 


Van Lieshout (2011) established for a general stationary point process 
the approximation J(r) — 1 ~ — p{K{r) — uJdT d } ) where ujd = \B(o, 1)| and 
r t —y UdT d is the AT-function for a stationary Poisson process. It is therefore 
not so surprising that often empirical J and A'-functions lead to the same 
practical interpretations. In particular, if for our LGCP c > 0, i.e., g > 1, 
then we have K(r) — 0Jd r<i > 0, and so we expect that J(r) < 1. Indeed Van 


Lieshout (2011) verihed this in the case of an LGCP with g > 1. This result 


immediately follows by the new expression (13) 


4.2 Laplace approximation 


Since Laplace’s pioneering work (see e.g. Stigler 1986), Laplace approxi¬ 


mations of complex integrals have gained much attention in probability and 
statistics, in particular when considering integrals involving Gaussian random 
fields 


^see e.g. 


Rue et al. 2009). This section discusses a Laplace approxima¬ 


tion of 1 — G(r); a Laplace approximation of 1 — F(r) can be obtained along 
similar lines. 

For A > 0, consider a grid of quadrature points, 


G{ A) = {(Aii, • • •, A i d ) | e Z}, 

and for v G G( A), let A^ = [iq — A/2, Vi + A/2[x ■ ■ ■ x [vd — A/2, Vd + A/2[ 
be the grid cell associated with v. Then for any non-negative Borel function 
i : —» M, we use the numerical quadrature approximation 


' B(o,r) 


exp{V (x)}£(x) dx 


y wJ(v)exp{Y(v)}, 

»S5(A)nB(o,r) 


(14) 
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where the quadrature weight w v = | fl B(o, r)|. 

Denote by /, M and £ the density, the mean vector and the covariance 


matrix of the normally distributed vector {F(u)}„ e g(A)n.B(o,r)- Then (12) and 
(|l4|) give 


1 - G(r) ~ / exp l - Y w v g(v) exp (y v ) l f(y) d y 

(«65(A)nB(o,r) 

= / exp{h(y)}dy 


(15) 


where y is the vector (y v )v&g(A)nB{o,r) of dimension m = #{£?(A) fl B(o,r)} 
and 

h(y) = ~ Y. «'„j(t>)exp(y„)-i( ! ,-M) T E- 1 (y-M)-llog{(2 W nS|}. 

v£g(A)nB(o,r) 

The gradient vector for h is 

Vh{y) = -d{y)-?T\y~M), (16) 

where d(y) = { w v g(v ) exp(^)}„ e g(A)ns(o,r ) 5 and minus one times the Hessian 
matrix for h is 

H(y) = D(y) + £"\ 

where D(y) is the diagonal matrix with entries {d(y)} v , vQ( A) fl B(o,r)}. 
Since H{y) is a positive definite matrix, h has a unique maximum at a point 
y, which can be found using Newton-Raphson iterations 


</' = yff> + II 1 {t/"Ui/' 1 }. 


(17) 


Therefore, the logarithm of the Laplace approximation of the right hand side 
in (15) (see e.g. Stiglcr, 1986) gives 


log{l - G{r)} « - w vg{y) exp (y v ) + ^{y- M) T d(y) 

v£Q(A)nB(o,r) 

-llog|D(i)E + /|. 


(18) 


where / is the m x m identity matrix. For the computation of £ 1 (y — M ) 
in (16) we solve LL T z = y — M where L is the Cholesky factor of £. In the 
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same way, considering the QR decomposition, Q(y)R(y) for y G M m , of the 
matrix D(y )E + /, the computation of H~'' 1 >{y^}'Vh{y^} in (17) is done 
by first solving Q{y^}R{y^}z = Vh{y^} and second by evaluating Sz. 
Finally, in Q, \D(y)E + 1 1 = \R(y)\. 


4.3 Numerical illustration 


To illustrate the Laplace approximations of the G and J-functions (Sec¬ 
tion 4.2) we consider three planar stationary LGCPs with intensity p = 50 


and spherical covariance function 


c(x) 



if ||x|| < a, 


0 


otherwise, 


with variance a 2 = 4 and scale parameters a = 0.1, 0.2, 0.3, respectively. We 
evaluate the approximations of G(r) and J(r) at r G 77, where 1Z is the set 
of 50 equispaced values between 0.01 and 0.25. For r G 77, we define the grid 
Q( A r ) with A r = 2 r/q, where q is a fixed integer. Such a choice implies that 
#{G( A r ) ft [—r, r] 2 } = g 2 , and so we have at least q 2 ir/4 quadrature points in 
B(o,r). For a given q, denote by G q , F q , and J q the corresponding Laplace 
approximations of G, F, and J, respectively. Figure [l] shows the resulting 
curves with q = 16. To see how far these Cox processes deviate from the 
Poisson case (which would correspond to a 2 = 0), we also plot the G -function 
in the Poisson case, namely 1 — G(r) = exp(— pur 2 ). To study the role of 
q, we report in Table [l] the maximal differences max r6 K |Gi 6 (r) — G q (r)\ and 
maXr eK \Ji e (r) — J q (r) \ for q — 4, 8 ,12. As expected, each difference decreases 
as q increases and is already very small when q — 12 (less than 4 x 1CF 3 except 
for the J-function and a = 1). This justifies our choice q = 16 in Figure [I] 
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G function 



I— 

a=0.3 

A - 

a=0.2 

- -O - 

o=0.1 

— 

Poisson 



'LA... 

- a.a.a 


0.00 0.05 0.10 0.15 0.20 0.25 


Figure 1: For three planar stationary LGCPs with intensity p = 50 and 
c given by a spherical covariance function, with variance a 2 = 4 and scale 
parameters a = 0.1, 0.2, 0.3, respectively, Laplace approximations of the G- 
function (left) and the J-function (right). 




q = 4 

q = 8 

q — 12 

a 

= 0.1, G 

59.9 

8.4 

2.1 


J 

505.9 

96.1 

20.5 

a 

= 0 . 2 , G 

14.3 

1.6 

0.5 


J 

109.0 

13.8 

3.5 

a 

= 0.3, G 

4.2 

0.5 

0.1 


J 

22.1 

3.1 

0.3 


Table 1: For the same three LGCPs as in Figure [lj maximal differences be¬ 
tween the Laplace approximations G q and G \ 6 . and between the Laplace 
approximations J q and J 16 , with q = 4, 8 ,12, that is, the evaluation of 
max rg 7 ^ |L/i 6 (r) — H q {r )| for H = G, J. Results are multiplied by 10 3 . 

The Laplace approximation of the G-function could also be derived using 
0- To check the agreement of the numerical approximations based on (JTTJ) 
and ( [12] ), respectively, Table [2] shows the maximal difference between the 
two approximations of first the G-function and second the J-function. In 
agreement with the theoretical developments, in both cases, the difference 
does not exceed 4 x 10 ~ 4 when q = 16. 
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q = 4 

q = 8 

q — 12 

q — 16 

a 

= 0.1, G 

3.8 

8.4 

4.4 

3.2 


J 

15 

8.5 

6.1 

3.9 

a 

= 0.2, G 

4.7 

3.1 

3.9 

1.5 


J 

4.7 

3.1 

2.1 

1.7 

a 

= 0.3, G 

0.1 

1.9 

1.4 

1.1 


J 

0.2 

2.0 

1.5 

1.3 


Table 2: For the same three LGCPs as in Figure [lj maximal differences 
between the Laplace approximations of the G -function based on (11) respec¬ 
tive (12), with q = 4,8,12,16, and similarly for the J-function. Results are 
multiplied by 10 4 . 


4.4 Scots pine saplings dataset 

The left panel in Figure [2] shows the locations of 126 Scots pine saplings in 
a 10 by 10 metre square. The dataset is included in the R package spatstat 


as f inpines, and it has previously been analyzed by Penttinen et al. (1992), 


Stoyan and Stoyan (1994), and Mpllcr et al. (1998). In the first two papers 


a Matern cluster process is fitted, using the K -function (or its equivalent L- 
function) and its nonparametric estimate both for parameter estimation and 
model checking, while the third paper considered an LGCP with exponen¬ 
tial covariance function and used the pair correlation function for parameter 
estimation and the F and G-fiuictions for model checking. Mpller et ah 


(1998) concluded that both models provide a reasonable fit although when 
also including a third-order functional summary (i.e., one based on Xj^) the 
LGCP model showed a better fit. Below we extend this analysis by using the 


J-function and the approximation established in Section 4.2 


We fitted both models by minimum contrast estimation (method kppm in 
spatstat) which compares a non-parametric estimate of the A'-function with 
its theoretical value. When approximating the J-function for the LGCP, we 
used q — 12; no improvements were noticed with higher values of q. The 
right panel in Figure [2] shows the theoretical J-functions for the two fitted 
models together with a non-parametric estimate of the J-function obtained 
from data, considering 50 equispaced distances (r-values) between 0 and 
0.9 meter (for the exact expression of the J-function for the Matern cluster 
process, see e.g. Mqllcr and Waagepetersen, 2004). Clearly, the fitted LGCP 
provides a better fit than the fitted Matern cluster process. Indeed, the 
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maximal difference between the non-parametric estimate and the theoretical 
J-function equals 0.43 for the Matern cluster model and 0.20 for the LGCP 
model. 




Figure 2: Left panel: Locations of 126 Scots pine saplings in a 10 by 10 
metre square. Right panel: Non-parametric estimate of the J-function (solid 
curve) and fitted J-functions for the Matern cluster process (dashed curve) 
and the LGCP with exponential covariance function (dotted curve). 


5 Concluding remarks 


We expect that our results for the reduced Palm distributions for an LGCP 
can be exploited further regarding third-order and higher order functional 
summaries (one such characteristic was briefly mentioned in Section 4.4), 
parameter estimation procedures, model checking, etc. This is discussed in 
some detail below. 

For likelihood based inference, suppressing in the notation any depen¬ 
dence on a parametric model and assuming a realization x of an LGCP is 
observed within a region W C K d of Lebesgue measure \W\ G (0, oo), the 
likelihood function is given by the density 


/(x) = E exp 


\W\ — / exp{y'(w)} dw 


'W 


Qexp{F(u)} 


ugx 
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with respect to the unit rate Poisson process on W. This density expression 
has no explicit form even for simple covariance function models for the under¬ 
lying Gaussian random held. For maximum likelihood estimation and pre¬ 
diction of the Gaussian random held, rather elaborate and time-consuming 
Markov chain Monte Carlo (MCMC) methods (Mpllcr and Waagepetersen, 
2004) may be used; alternatively, a Bayesian approach based on integrated 


nested Laplace approximations or MCMC methods may be used (Rue et al. 


2009 Taylor and Diggle, 2014). Using the results of this paper, we have the 


following new expression for the density: 


/(x) = p (n) (x i,... ,x n )E fexp 


\W\ 


exp{y; 


XI, 


(w)} dw 


'w 


Here the expression for the n-th order intensity is explicit, but it remains 
to investigate if the expectation, which apparently has a simpler expression, 
may be easier to approximate. 

For model checking, when considering non-parametric estimates of func¬ 
tional summaries together with simulated confidence bounds under a htted 


LGCP model (such as extreme rank envelopes, see Baddeley et al. 2015 


Myllymaki et ah, 2016), it could be pertinent to include the theoretical ex¬ 


pressions of the functional summaries for LGCPs obtained in the present 
paper. 

Finally, recall that for any point process, the pair correlation function 
(when it exists) is invariant under independent thinning. Could this property 
be exploited in connection to LGCPs where we know how the pair correlation 
function is related to those of the reduced Palm distributions? 
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A Proofs 


Proof of Lemma [TJ- By conditioning on A, (|TJ) becomes 

7 L 

^ h(x 1 ,...,x n ,X\{x 1 ,...,x n }) 


EE 


A 


EE I J ' J K x u ■ ■ • , x n, X) 11 A(xj) dx 

J ■■■ J E ... ,z n ,X) JjA(xi)J 


1 ■ ■ ■ dx r 


A 


d^i• ■ • dx r 


(19) 


( 20 ) 


Here, in (19) we use that X given A is a Poisson process and apply the 
extended Slivnyak-Mecke theorem (Mpller and Waagepetersen 2004), and 
in (20) we use Fubini’s theorem. Combining g and (20), we deduce g. 
Finally, g follows from g with h(x i,..., x n , x) = l(x D K = 0). 

Proof of Theorem 0 By g and <@-0. we just have to show that for any 
compact K C S and pairwise distinct points x±,... ,x n G S, 


Eexp 
= Eexp 


exp 


' K 


u) + n Xl ,...,xM \ du 


£ Y(xi) - c(xi, Xj)/ 2 - / exp <j n(u) + Y(u) \ d u 

i,j =1 


i=l 


' K 


This will follow by verifying that the distribution of {Y(-u) + ^’ 1 =1 c(u, (Cj)} ue ,s 
is absolutely continuous with respect to the distribution of Y = {Y(rt)} ue s, 
with density exp jX^Li v( x i) ~~ Y^j=i c ( x i> x j)/^f when y is a realization 

of Y. Since the distribution of a random held is determined by its finite 
dimensional distributions, we just need to verify the agreement of the char¬ 
acteristic functions of the probability measures Q i and Q 2 given by 


Qi(B) = P 


f ~ n 
Y(U) + Y, 
^ 1=1 


c(u,Xi) 


e B 


u£U 


and 


Q 2 (£)=E 1 {Y(u)} u£U G B 


exp IX Y(xf) y ^ c^XiiXj )/2 
i,j = 1 


i =1 
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for any Borel set B C M n+m , any pairwise distinct locations u \,..., u rn e 
R d \ {x\,... ,x n }, with m > 0 and U = {x \,..., x n ,ui,... ,u m }. Let E = 
{c(u, v)} UjVe u denote the covariance matrix of {Y'(u)} u£ u, and let c Xl) ... )Xn = 
(SiLi C ( M ’ x i)}u&u — Ele, where e consists of n l’s followed by m 0’s. For 
t G M n+m , the characteristic function of Q 2 is (with i 2 = —1) 


Eexp \ i{Y(u)}Z eI t + {Y(u)}l eI e + e T Ee/2 
= exp (e T Ee/2) Eexp i{Y — ie) 

= exp (e T Ee/2 + ie T St — t T Et/2 — e T Ee/2) 
= exp (ie T St — t T Et/2) . 


The last expression is the characteristic function of Q\ which concludes the 
proof. For the second last equality in the above derivation we considered 
{Y'(u)} u£ u as a complex Gaussian vector and used the expression for its 
characteristic function. 
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